Binary black hole mergers in gaseous disks: Simulations in general relativity 
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Simultaneous gravitational and electromagnetic wave observations of merging black hole binaries 
(BHBHs) can provide unique opportunities to study gravitation physics, accretion and cosmology. 
Here we perform fully general-relativistic, hydrodynamic simulations of equal-mass, nonspinning 
BHBHs coalescing in a circumbinary disk. We evolve the metric using the Baumgarte-Shapiro- 
Shibata-Nakamura (BSSN) formulation of Einstein's field equations with standard moving puncture 
gauge conditions. We handle the hydrodynamics via a high-resolution shock-capturing scheme. 
These initial simulations are exploratory in nature and simplified accordingly. We track the inspiral 
starting from a binary separation of lOM, where M is the total binary mass. We take the disks 
to have an inner radius at Rin « 15A4" to account for the hollow created by the binary torques. 
Our disks extend to i? « 65Af and have an initial scale height of H/ R « 0.03 — 0.11. The gas is 
governed by a F-law equation of state, with F equal to 5/3, 4/3, and 1.1. Disks are allowed to relax 
in the "early inspiral" epoch to provide quasistationary realistic initial data. We then evolve the 
spacetime metric and matter during the "late inspiral and merger" epochs. The later simulations are 
designed to track BHBH inspiral following disk-binary decoupling, through merger and ringdown, 
terminating before viscosity has time to fill the hollow about the black hole remnant. We compute 
the gas flow and accretion rate and estimate the electromagnetic luminosity due to bremsstrahlung 
and synchrotron emission as a perturbation for optically thin disks. The synchrotron component of 
the luminosity peaks in the infrared band and should be detectable by WFIRST and possibly the 
LSST for a 10^ Mq binary embedded in a disk with a density n ~ lO'^^ cm"'' at 2 = 1, beginning 
with a maximum value of i ~ 10'*^nf2Af|erg s~^ at decoupling, and decreasing steadily over a 
timescale of ~ 100 Mg hours to a value of L ^ 10*^ 71^2 Af| erg s~^ at merger. 

PACS numbers: 04.25.D-, 04.25. dg, 47.75. +f 



I. INTRODUCTION 

All bulge galaxies (including the Milky Way) are 
believed to contain a central supermassive black hole 
(SMBH) with a mass M between 10"* M© and IO^Mq [J- 
y] . It is also believed that galaxy mergers commonly lead 
to the formation of a massive black hole binary (BHBH) 
system in the merged remnants [1, Q. In the stan- 
dard picture, the BHBH binary separation decreases, first 
through dynamical friction due to distant stellar encoun- 
ters, then through gravitational slingshot interactions in 
which nearby stars are ejected at velocities comparable to 
the binary's orbital velocity, and finally through the emis- 
sion of gravitational radiation, leading to coalescence of 
the binary [Q]. These low-frequency gravitational waves 
will be detectable by LISA and will contain a wealth of in- 
formation about the inspiral. The gaseous accretion flow 
that forms around the binary can be a source of electro- 
magnetic radiation as well. Following the detection of 
gravitational waves from a BHBH merger, electromag- 
netic "afterglow" radiation can provide confirmation of 
the coalescence The timescale during which de- 

tectable "afterglow" radiation rises to its maximum is 
governed by viscous diffusion of gas close to the remnant 
and ranges from several years to hundreds of decades in 
the case of supermassive BHBH systems. 
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There also exists the possibility of detecting electro- 
magnetic "precursor" radiation prior to the merger and 
before the maximum gravitational wave emission |l^.[l5j. 
If the merger takes place in a hot gas cloud in which 
the distant gas is nearly homogeneous and either at rest 
with respect to the binary ("binary Bondi" accretion) 
or moving ( "binary Bondi- Hoyle-Lyttleton" accretion), 
then Farris et al. [la] (hereafter Paper I) have shown that 
the luminosity will peak at the end of the binary inspiral 
phase immediately prior to the final plunge. At this stage 
shock heating of the gas and turbulent magnetic field 
amplification are strongest. The peak luminosity lasts 
for 5t ^ Mq hours prior to merger and then plummets 
sharply following the coalescence. Here Mg is the binary 
mass in units of lO^M©. If, instead, the accretion takes 
place via a geometrically-thin, optically-thick Keplerian 
disk around the binary ("binary Shakura-Sunyaev" ac- 
cretion), there may be a late-time precursor brightening 
from tidal and viscous (or turbulent magnetic) dissipa- 
tion in the inner disk. This radiation peaks on a timescale 
of 6t ~ O.lMg days prior to merger and remains high af- 
terwards [Til . Each of these scenarios raises the exciting 
possibility of a simultaneous detection of electromagnetic 
and gravitational waves from a BHBH merger. 

This picture is loosely supported by a number of ob- 
served AGNs that may be harboring BHBH binaries. 
Very-long baseline interferometry (VLBI) observations of 
the elliptical galaxy 0402-1-379 have discovered two radio 
sources at a projected separation of only 7 pc. The ex- 
istence of jets, as well as variability associated with BH 
activity, indicate that the system may be a BHBH bi- 
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nary |17H19l |. Another candidate is OJ 287, a BL Lac 
object whose hght curve shows variabihty with a period 
of ~ 12 yr. It is beheved that this may be a massive 
BHBH binary around which the smaUer BH orbits with 
a period of 12 yr, penetrating the disk of t he p rimary 
and giving rise to the observed variabihty [20l - l22l | . It has 
been proposed that the quasar SDSS 092712.65+294344 
may be either a binary system [23l [23 |. or a recoil- 
ing BH that is the product of a binary merger p5| . 
Such suggestions are supported by a systematic shift of 
2650 km s^^ in the emission lines. Another candidate is 
quasar SDSS J153636.22-I-044127.0, in which two broad- 
line emission systems are observed, separated in velocity 
by 3500 km s^^ . This observation has been interpreted as 
a BHBH binary system in which each object has its own 
emission system [2g|. Recently, the first triple AGN sys- 
tem, SDSS J1027-hl749, has been discovered [ijj. This 
galaxy contains three emission-line nuclei corresponding 
to a supermassive black hole triple with kpc-scale sepa- 
rations. 

Information from a simultaneous detection of elec- 
tromagnetic and gravitational waves may be useful for 
studying fundamental aspects of gravitational physics. 
For example, in some modified gravity scenarios, the 
propagation velocity for gravitons may differ from that 
of photons [1^, [2^ . Additionally, the measurement of the 
luminosity distance from the gravitational wave signal at 
an accuracy of 1% — 10%, coupled with the redshift infor- 
mation from the electromagnetic detection, could serve 
as a cosmological "standard siren" of unprecedented ac- 
curacy (better than ~ 1%) [s^l- Such detections may also 
combine accurate measurements of BH spins and masses 
obtained from gravitational wave signals with electro- 
magnetic observations to probe BH accretion physics in 
great detail 3l|. It has even been proposed that simul- 
taneous detections of electromagnetic and gravitational 
waves r nay provide a means of witnessing the birth of a 
quasar [32| . 

Several mechanisms for electromagnetic emission from 
accretion disks around merging BHBH binaries have been 
proposed. In one scenario, the inner edge of the accre- 
tion disk is identified as the radius at which the viscous 
torque on the gas balances the gravitational torque from 
binary. This radius is between 1.5 and 2 times the or- 
bital separation [33|-(36{ and encompasses a hollow region 
in the disk. Late in the inspiral the BHBH binary decou- 
ples from the disk and coalesces. For a binary of mass 
M K. 1O^M0 accreting at 10% of the Eddington rate, the 
subsequent evolution of this disk, which is optically thick, 
gives rise to a source that initially peaks in the UVh&nd 
and then hardens to extreme ultraviolet and soft x-ray 
emission at late times 0, [H, [l^ ■ Additionally, there is a 
sudden change in the mass of the binary during the final 
stage of the merger, as gravitational waves carry away a 
few percent of the mass. The abrupt change in the poten- 
tial creates waves in the disk which may grow into shocks 
and increase the luminosity of the disk in a unique way 
giving rise to a detectable prompt x-ray sig- 



nal. Another possibility is that the merged BH remnant 
may experience a recoil velocity which c an, in principle, 
be as high as several thousand km s"-"^ [S^, although it 
is likely to be much lower (< 200 km/s) in most galaxy 
mergers [s^. This recoiling BH may "shake" or pene- 
trate the disk, creating shocks which could give rise to a 
transient signal. 

Various methods have been used to model plausible 
sources of electromagnetic emission from BH mergers. 
In one approach, the dynamics of the inspiral is ignored, 
focusing on the effect of BH kicks and/or BH mass loss 
on the hydrodynamical flow [sl-flll l37l [jol - fi^ . In another 
approach, the behavior of the gas is modeled by following 
the motion of collisionless "particle tracers" on geodesies 



44| . Other approaches involve vacuum and/or force- free 



calculations to investigate the role that magnetic fields 
may play in producing detectable electromagnetic emis- 
sion when the density near the binary at merger is very 
low [4^ li^. Only recently have fully relativistic, hy- 
drodynamical simulations of BHBH binary inspiral and 
merger in a gaseous environment been performed [l6l . IZtI - 

Si. 

In this paper we study BHBH binary mergers in the 
presence of a circumbinary gaseous disk. Modeling such 
systems requires fully general-relativistic dynamical sim- 
ulations. The development of stable algorithms to inte- 
grate Einstein's field equations of general relativity nu- 
merically in 3 -|- 1 dimensions, such as the Baumgarte- 
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Shapiro-Shibata-Nakamura (BSSN) formalism ^ ^ 

and the generalized harmonic approach [52h5^ , together 
with the identification of suitable gauge conditions, has 
enabled several pioneering simulations that demonstrate 
how to track the late inspiral, plunge and merger of a 
BHBH binary in vacuum [55| - |57| . More refined follow- 
up simulations of these strong-field, late phases, joined 
onto analytic, pos t- Newtonian calculations of the early 
inspiral epoch (58| . are now capable of producing accu- 
rate gravitational waveforms for merging BHBH binaries 
with companions spanning a range of mass ratios and 
spins (see e.g. [5^ and references therein). 

With the problem of gravitational wave emission from 
a vacuum BHBH binary inspiral well in hand, it is now 
important to turn to the problem of electromagnetic 
emission from BHBH binary coalescence in an astrophys- 
ically realistic gaseous environment. 

In Paper I, we considered hot accretion flows in which 
the gas is near the galaxy virial temperature and the spe- 
cific angular momentum of the gas L is less than that of 
a circular orbit near the horizon, L ^ Mc. Such fiows 
can be well described by the spherical Bondi, or Bondi- 
Hoyle-Lyttleton accretion model. In this paper, we con- 
sider flows in which the angular momentum of the gas 
cannot be neglected, and the accretion flow is disklike. 

Disk accretion onto a BHBH binary has been stud- 
ied previously in the Newtonian, geometrically thin-disk 
limit, both anal ytic ally [l^ [l^ [isl . |60-62] and numeri- 
cally [1, [T^, [ill, l36j | . We extend this work by perform- 
ing fully relativistic hydrodynamcal simulations in 3-f 1 



3 



dimensions. Our treatment here is quite preliminary 
and meant to introduce the computational framework 
for more detailed and realistic simulations that we are 
preparing. In this paper we restrict our attention to a 
circumbinary disk residing in the orbital plane of two 
nonspinning, equal-mass, binary black holes. The black 
holes are initially in a quasistationary, nearly circular 
orbit and represent a solution to the conformal thin- 
sandwich (CTS) initial-value equations (see, e.g. j63l - [65j 
and references therein) . The mass of the disk is assumed 
to be small in comparison to the total black hole mass. 
We explore the response of the disk to the binary on 
a dynamical timescale and thus ignore the secular mo- 
tion of gas in the disk due to viscosity or turbulent mag- 
netic fields. We treat the gas as a perfect fluid described 
by a F-law equation of state (EOS) and handle shocks 
by high-resolution, shock-capturing (HRSC) techniques. 
We study the response of the disk to tidal torques dur- 
ing the early and late inspiral phases, as well as during 
the merger and post-merger epochs. The inspiral and 
merger are followed by solving the BSSN field equations 
[so, 51] with moving puncture gauge conditions 56, 57]. 
We are particularly interested in the evolution of the hol- 
low region in the disk about the binary [33l - l36j and the 
extent to which gas penetrates the hollow and accretes 
onto the black holes. We also estimate, as a perturbation, 
the electromagnetic emission from the disk that charac- 
terizes the inspiral and merger epochs. Our treatment 
is appropriate for describing the epoch following disk- 
binary "decoupling," when the BHBH inspiral timescale 
is much shorter than the viscous timescale in the disk, 
whereby viscosity-induced inflow can be neglected. Our 
analysis remains valid throughout the binary merger and 
ringdown phases, but is no longer adequate to describe 
the late-time evolution when viscosity serves to drive gas 
into the hollow and accrete onto the merged remnant 
0,111. We briefly compare our results with another, 
recently published, general relativistic study [i^ that 
treats a similar scenario, but employs different methods 
(e.g. different initial data and luminosity estimates) and 
addresses different issues. 



The structure of the paper is as follows. In Sec. [IT] 
we summarize the computational challenge posed by the 
large dynamic range associated with this problem, and we 
describe our approach for overcoming this challenge. In 
Sees. Illll and llVl we briefly outline the basic gravitational 
field and matter evolution equations and their specific 
implementations in our numerical scheme. Here we also 
provide an overview of our initial data, gauge conditions, 
and diagnostics. In Sec. IIVDI we review code tests that 
were performed to validate our numerical scheme. In 
Sec.|V]we describe the results of our binary BHBH merger 
simulations. In Sec. IVII we summarize our findings, and 
briefly compare with previous simulations. Henceforth 
we adopt geometrized units and set G = c = 1. 



II. COMPUTATIONAL CHALLENGE 

Simulating realistic accretion flows is extremely chal- 
lenging due to the enormous dynamic range character- 
izing the time and length scales in the problem. One 
length scale is set by the initial, total Arnowitt-Deser- 
Misner (ADM) mass of the binary system, M . We ne- 
glect the mass of the disk and assume that Mdisk ^ M . 
The ADM mass sets the length scale at which rela- 
tivistic effects become significant. Another important 
length scale is the binary orbital separation a. Asso- 
ciated with the orbital separation is the orbital period, 
torb = '2i^/^bin ~ 27r(a/Af)"^/^, where ^lun is the binary 
orbital angular velocity. 

Torques from the binary have the effect of driving mat- 
ter in the vicinity of the BHBH orbit outward, creating a 
hollow cavity inside the inner edge of the accretion disk 
at Rin ^ 1.5 — 2a [33l - l36j . This radius is determined by 
the balance between viscous stresses in the disk and tidal 
torques from the binary [U |3^ , provided that the vis- 
cous timescale t^s ~ (2/3)i?f„/i^ is shorter than the grav- 
itational inspiral timescale tmerge ~ (5/16)a'*/M'^. Here 
ly is the shear viscosity, and we assume an equal-mass 
binary. As the orbital separation decreases, the binary 
eventually enters an epoch at which tmerge < tms- At this 
point, the binary decouples from the disk (tI. [Til [63. l66j . 
If one assumes an a-disk with a viscosity law ^{R) = 
{2/3)aPgas/{p^K), where p is the gas density and Pgas 
is the gas pressure, then the decoupling radius a^ is given 
by UMM 

^ « 126a:r/'°5-49/200A^/i°Af6'/''(<5-iC)^'/'°ea2'^'°" , 

where a — 0.1a_i, S = 0.l5_i, S = 37rI](ad)i^(ad)/-^'^Edd 
and 9 = 0.26'o.2- Here MEdd = A-KMrUp/ {-qar) is the Ed- 
dington accretion rate, ctt is the Thomson cross section 
for electron scattering, 77 is the radiative efficiency, 9 <1 
is a porosity correction factor applied to the scattering- 
dominated optical depth [67J , and 5 roughly accounts for 
the shortening of the viscous timescale at the disk edge 
where the surface density S is very steep [68|. 

Another important length scale is the characteristic 
size of the disk, Rdisk , which we define here as the radius 
at which the gas pressure is maximum, Rdisk = R{Pmax)- 
In general, Rdisk depends on the details of the temper- 
ature and angular momentum profile in the disk, and 
is highly dependent on the particular choice of disk 
model. Associated with Rdisk is the orbital time scale 
idisk which we define as the Keplerian orbital period 
td^sk = 2^:{Rd^sk/Mf/^. 

If we assume the size of the entire disk is 
several Rdisk 3> Rin and use the estimate of Od given in 
Eq. ([1]) , we find that a simulation of the full inspiral from 
decoupling to merger would require us to resolve length 
scales from ~ M to ^ 10"^M. More challenging, we must 
resolve timescales from M to several tmerge ^ l(fM. 
Unfortunately, the latter is beyond the capability of cur- 
rent numerical codes. In order to circumvent this is- 
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sue, we consider a disk with relatively small values of 
ttd = lOM, ~ 15M, and Rd^sk ~ 35M. With these 
choices, the important time scales become torb = 225M, 
tdisk = 1300M, and tmerge — 1250A'/. Givcn the wide 
range of gaseous environments in galactic cores, such pa- 
rameters are not implausible, and we expect that qualita- 
tive features of our results can be extended to disks with 
larger values of ad, Rin, and Rdisk- Our choice allows 
us to study the full evolution of the system from decou- 
pling to binary merger. We focus on the post-decoupling 
phase through merger, ringdown, and disk equilibration, 
but prior to disk inflow on viscous time scales. Accord- 
ingly, our perfect-fluid approximation will closely mimic 
a realistic flow during these epochs, as the viscous time 
scale (which may originate from MHD turbulence) is long 
compared to the length of our simulations. 



III. BASIC EQUATIONS 

Throughout this paper, we use Latin indices to denote 
spatial components (1-3) and Greek indices to denote 
spacetimc components (0-3). 



A. Early inspiral epoch 

We define the "early inspiral epoch" as the phase of 
the binary inspiral prior to decoupling. Throughout this 
phase, the inspiral time scale is much longer than the or- 
bital time scale. This fact can be exploited by neglecting 
the change in binary separation and employing a metric 
that is quasistationary in the rotating frame of the binary. 
This simplification provides an accurate solution for the 
spacetime without the computational expense of a full 
evolution of Einstein's field equations. We evolve the full 
relativistic hydrodynamics equations in this background 
metric over ~ 5tdisk to enable the disk to relax to a qua- 
sistationary state. This technique thus provides astro- 
physically realistic initial data with which to begin evo- 
lution of the late inspiral and merger epochs fSec lIIIBjl . 

In order to use this method, we must choose a coor- 
dinate system in which the metric explicitly refiects the 
symmetry of the spacetime. This symmetry, describing 
a spacetime that is quasistationary in a frame that ro- 
tates with the orbital frequency of the binary ft, can be 
constructed by employing a helical Killing vector, 

^ = dt + nd^. (2) 

For a spacetime admitting such a Killing vector, we have 
^eg^iiy = , (3) 

where C is the Lie derviative, and is the spacetime 
metric. 

Provided we are working in an appropriate coordinate 
system (i.e. one employing Killing coordinates t and 0), 



we may express the metric at any point in spacetime in 
terms of the metric on an initial t = slice according to 

g^,{t,r,0,(f>)=gf,,iO,r,e,c^-nt) (4) 

where the symbol = denotes that the equality holds only 
in a particular coordinate system. One can easily verify 
that the above equation satisfies Killing's equation 

We note that Eq. ^ is written in spherical polar co- 
ordinates, i.e. {a;"} = {t,r,6,(j)}. However, Cartesian 
coordinates are more suitable for work in 3D, as coor- 
dinate singularities at r = and on the polar axis are 
avoided. We therefore transform the spherical compo- 
nents of g^^ back to the Cartesian components using the 
usual tensor transformation formula. 

BHBH evolution employing standard puncture initial 
data and moving puncture gauge conditions does not re- 
sult in a metric that satisfies Eq Q (Puncture initial 
data does not implement a helical Killing vector). By 
contrast, BHBH CTS initial data (see, e.g. [64|) are 
specifically constructed to satisfy this equation: CTS ini- 
tial data impose the condition that the spacetime in the 
rotating frame is stationary (see [65] for discussion and 
references). This condition is valid, approximately, when- 
ever the binary companions are sufficiently well separated 
that the inspiral time scale is much longer than the or- 
bital time scale. In this quasistationary early inspiral 
regime we can employ CTS initial data and CTS lapse 
and shift functions to evolve the metric via a simple co- 
ordinate rotation in lieu of integrating the Einstein field 
equations. We can then evolve the disk by integrating the 
hydrodynamic equations for the fluid in this background 
spacetime. 

CTS initial data contains excised interiors. We fol- 
low the technique outlined in [69i] and fill the excised re- 
gion inside the BH interiors with smoothly extrapolated 
"junk" (i.e., constraint-violating) data. This treatment 
is valid because the interior regions are causally discon- 
nected from the exterior spacetime. 

B. Late inspiral and merger epochs 

We evolve both the metric and hydrodynamic equa- 
tions during the late inspiral and merger epochs. Our 
basic gravitational field and relativistic hydrodynamics 
equations are discussed in [sil . Fzoj , where their numerical 
implementation is described and detailed code tests are 
summarized. Here, we briefly review these equations and 
their implementation. 

We write the spacetime metric in the standard 3-1-1 
form, 

ds^ = -a^dt^ + 7y- {dx' + I3'dt){dx^ + P^dt). (5) 

where a, /3*, and 7y are the lapse, shift, and spatial met- 
ric, respectively. The extrinsic curvature Kij is given by 

{dt - Cp)^^j = ~2aK,j, (6) 
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where £^ is the Lie derivative with respect to B\ We 
evolve 7ij- and Kij using the BSSN formulation |50l. [5l|. 
The fundamental variables for BSSN evolution are 



= 


— ln[det(7y)] , 


(7) 


7»i = 




(8) 


K = 




(9) 


Aij = 




(10) 


r = 




(11) 



The evolution and constraint equations for these fields 
are summarized in [s^, IMI- We assume in this paper 
that the mass of the gas is negligible compared to the 
mass of the BHs, thus we do not include matter source 
terms in our metric evolution equations. 

Adding Kreiss-Oliger dissipation to the BSSN evolu- 
tion equations outside the BH can reduce high-frequency 
numerical noise associated with adaptive mesh refine- 
ment (AMR) refinement interfaces [s^ • We use this tech- 
nique here and have found it useful in reducing Hamilto- 
nian and momentum constraint violations. 

We adopt the standard moving puncture gauge con- 
ditions: an advective "l-|-log" slicing condition for the 
lapse and a "Gamma-freezing" condition for the shift 
[711 . Thus we have 

doa = -^2aK, (12) 

doP' = (3/4)i?^ (13) 

doB' = dot' - r^B' , (14) 

where do = dt ~ P''dj. The rj parameter is set to 0.5/Af 
in all simulations. 

C. Evolution of hydrodynamic fields 

The fundamental matter variables are the rest-mass 
density po, specific internal energy e, pressure P, and 
four- velocity . The stress-energy tensor for an ideal 
gas is given by 

where h = 1 + e + P/po is the specific enthalpy. We 
evolve the "conservative" variables p,, Si, and f. They 
are defined as 

P* = -^/iPon^iU^ (15) 
~S^ = ^/iT^.n^-i", (16) 
f = ^Tf^un^n" - . (17) 

Here = (o!^^, — a"^/?') is the timelike unit vector nor- 
mal to the t — constant time slices. Evolution equations 
are given by Eqs. (34), (36), and (38) of 

dtp*+dj{p.v^)^Q, (18) 
dtS, + d,{a^Ti,) = ^a^T"f'^^go.0, (19) 
dtf + 9,(«2V7r°' - P*v') = s , (20) 



where 7 = det(7y) = e^^*^ and = /u^ is the fluid's 
3-velocity. The energy source term s is given by 

-I^T°°I3' + T°')^^a] . (21) 



D. Equation of state 

To complete the system of equations, we must specify 
an EOS. While our code can handle any EOS of the form 
P = P{po,e), we adopt a F-law EOS, 

P=(r-l)poe. (22) 

We perform simulations with F — 4/3, 5/3, and 1.1. By 
varying F we effectively examine gas flow under a full 
range of conditions. We choose F = 5/3 as our canonical 
case. The choice of F = 1.1 approximates an isother- 
mal gas (we have chosen F = 1.1 rather than F = 1 in 
order to retain the F-law form of the EOS while still ap- 
proximating isothermality). At t = 0, we take the EOS 
to be isentropic, with P — KpQ, where K = constant. 
Throughout this paper, we define temperature by 

P = 2nkT , (23) 

appropriate for pure ionized hydrogen. 



IV. NUMERICAL METHODS 
A. Disk initial data 

For our disk initial data, we use the equilibrium solu- 
tion for a stationary disk around a single Kerr BH derived 
by Chakrabarti et al. [tII and summarized in We 
take this disk as inital data for a binary BHBH, placing 
the inner boundary of the disk well outside the BHBH 
orbital radius. Though no longer stationary, the initial 
disk settles down to quasistationary equilibrium as the 
binary rotates, apart from low amplitude spiral density 
waves induced by the time- varying tidal torque. For com- 
pleteness, we provide a brief summary of the construction 
of disk initial data in Appendix Rl 

For our fiducial equation of state, F = 5/3, the re- 
sulting outer disk radius is Rout ~ 65M and the disk 
scale height at Rdisk is H/R = 0.11 (see TableUfor more 
details). Here H is defined as the height above the equa- 
torial plane where the pressure falls to 1/e its value on 
the equatorial plane at the radius of maximum pressure. 
For binary BHs, the disk is approximately stationary if 
Rin ^ a. Initially, we take Rin/a = 1.5. We find that 
the disk relaxes to a near quasistationary state after a 
time ~ 'it disk- 
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B. Evolution of metric and matter 

We evolve the BSSN field equations with fourth-order 
accurate, centered, finite-difference stencils, except on 
shift advection terms, where we use fourth-order accurate 
upwind stencils. We apply Sommerfeld outgoing wave 
boundary conditions to all BSSN fields. Our code is em- 
bedded in the Cactus parallelization framework [tJ , and 
our fourth-order Runge-Kutta time-stepping is managed 
by the MoL (Method of Lines) thorn, with a Courant- 
Friedrichs-Lewy factor set to 0.5 in all BHBH simula- 
tions. We use the Carpet [t^ infrastructure to implement 
the moving-box AMR. In all AMR simulations presented 
here, we use second-order temporal prolongation, cou- 
pled with fifth-order spatial prolongation. The apparent 
horizon of the BH is computed with the AHFinderDirect 
Cactus thorn (76| . 

We write the general-relativistic hydrodynamics equa- 
tions in conservative form. They are evolved by 
an HRSC technique [T^I that employs the piecewise 
parabolic (PPM) reconstruction scheme [73] coupled to 
the Harten, Lax, and van Leer (HLL) approximate Rie- 
mann solver [78|. The adopted hydrodynamic scheme 
is second-order accurate for smooth flows, and first- 
order accurate when discontinuities (e.g. shocks) arise. 
Throughout the evolution, we impose limits on the pres- 
sure to prevent spurious heating and negative values of 
the internal energy e. Specifically, we require Pmin ^ 
P < ^max inside the horizon, where Pmax = lO-ftTpg and 
-Pmin = Kpq/2. Whenever P exceeds Pmax or drops be- 
low Pmin, we reset P to Pmax or Pniin, respectively. We 
check that this fix is applied only inside the apparent 
horizon, which is causally disconnected from the rest of 
the grid. 

At each timestep, we need to recover the "primitive 
variables" po, P, and from the "conservative" vari- 
ables p*, T, and Si. We perform the inversion as speci- 
fied in Eqs. (57)-(62) of [701, but with a slightly modified 
analytic quartic solver from the GNU Scientific Library 
that outputs only the real roots. We use the same tech- 
nique as in [ll] to ensure that the values of Si and f yield 
physically valid primitive variables, except we reset f to 
10- '^'^Tg^inax (whcrc TQ^max IS the maximum value of t ini- 
tially) when either Si or f is unph ysic al (i.e., violate one 
of the inequalities (34) or (35) in [T^). The restrictions 
usually apply only to the region near the puncture inside 
the horizon. 

For each of our calculations, we set our outer bound- 
ary at 128Af and use 8 AMR refinement levels. The 
maximum resolution near each BH is Sx/M = 0.03125. 
For our single BH test calculations, we place our outer 
boundary at 128A'/ and use 6 AMR refinement levels. 
For these cases, the highest resolution near the BH is 
5x/M = 0.0625. 

We model the emission of electromagnetic radiation by 
treating this radiation loss as a perturbation, and neglect 
its influence on the hydrodynamic flow, as well as any 
deviation from adiabaticity that it induces. 



C. Diagnostics 

1. Surface density 

In order to track the global evolution of disk structure 
and compare with other disk calculations, it is useful to 
define the surface density E. Following [80], we define 

E(P, (/))=/ Pou*V^dz , (24) 

Jz>0 

where R = y^x^ + {R will always be the cylindrical 
radius in this paper, while r will always be the spherical 
polar radius) . We also report the angle-averaged surface 
density (I](P)) where 

(E) = — / Ed^ . (25) 

2. Flux diagnostics 

To derive meaningful fiux diagnostics we must iden- 
tify the conserved currents. Details of this derivation are 
given in Appendix A of (l6j . To summarize, consider a 3D 
region which lies between two world tubes F and L on 
a timeslice i=const. Let F be defined by h(t, x, y, z) — 0, 
and L be defined by l{t,x,y,z) = 0. In |16|, this region 
is depicted as the lower shaded region in Fig. 24. For 
the purposes of this paper, we let F be the world tube 
defined by the apparent horizon(s), and L be the world 
tube defined by a sphere of constant coordinate radius 
centered at the origin. Of course the surfaces could be 
chosen to take on other shapes as well. 

3. Conserved quantities 

Now consider a conserved current, j'^ which satisfies 
y^j^ = . (26) 
Then it can be shown that (see e.g. Appendix A of [T^]) 

q=^^-J'F+J'L . (27) 

where 

q{t) = [ fd^n^ (28) 

= / fV^ d?x , (29) 

J^p = - ( y^fdadb , (30) 
Jf 

^L = - j y^fdadb . (31) 
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Here g' is the determinant of the metric in the (t, h, a, b) 
or (t,l,a,b) coordinate systems. In the above example, 
is the flux of q across the horizon(s), while J^l is the 
flux of q across the outer sphere. 



4- Freedom in coordinate choice 

These fluxes are independent of any changes in coordi- 
nates that leave the slicing intact. Equivalently, we may 
alter the shift without affecting the flux, but the lapse 
must be kept the same. We can rewrite these fluxes in 
any other coordinate system (t, x, y, z) which preserves 
the same slicing. While a and b can be any two coor- 
dinates on the surface, we label them here as 9 and cf) 
for convenience, as this is done in our actual numerical 
calculations. 



J-p = I \f--g del 
If 



d{x,y,z) 



-g det 



d{h,e,^) 

d{x,y,z) 



8(1,0,^) 



fd^,h ded(t) , (32) 
fd^,l d9d(t) . (33) 



5. Rest-mass conservation 

Rest-mass conservation, Vp(pow^) = 0, corresponds to 
— Pqu^^. If we now define 



T 



Mq = I ^/—gpQu d x = / p^,d X 
d{x,y, z) 



(M) 



(M) 



^—g det 



-g det 



d{x,y,z) 



d{l,6,4>) 



(34) 

Pou^'h^^ d0d(j) ,(35) 
pouf^l^, dOdcf) , (36) 



then we may integrate Eq. (j27p in time to derive the 
following rest-mass conservation law 

(^Mo {t)+ dt (F^i'^ - J-f ^) ) ^ /Mo,, = 1 , (37) 

where Mg^i is the rest mass between the horizons and 
the surface L at i = 0. Equivalently, we can define the 
rest-mass accretion rate 



(38) 



6. E — Q.J conservation 



We employ the helical Killing vector defined in Eq. ([2]) 
to construct another conserved current. 



(39) 



We now define the following quantities 

Eit) = - [ d'x , 



Jit) = / V^T'^ d^x 



(E) _ 



{E) 



-niJ) = _ 



~g det 



-g det 



d{x,y,z) 



d{h,0,^) 
d{x,y,z) 



0(1,0, (b) 



T^'tK^ d0d(j) 



(40) 
(41) 

(42) 
(43) 



(J) 



-g det 



-g det 



d{x,y,z) 



d{h,0,(j)) 
d{x,y,z) 



d{l,0,cf) 



T%h^^ d0dc^ , (44) 
T"^?,^ d0d^ , (45) 



and we see that Eqs. \21\ and (p9l) give 



(46) 



for spacetinies possessing a helical Killing vector. Again, 
we may integrate in time to derive another conservation 
law 



E{t) - nj{t) 



dt ( ^L^) - nj^' 



dtlT 



7(E) 



E{t,) - nj{t,) 



(47) 



7. Spiral density wave diagnostics 



In our simulations, we use our CTS metric and Eq. (|4]) 
to ensure that our disk models exhibit quasistationary 
behavior before we begin the binary inspiral. Such qua- 
sistationary configurations are interesting in their own 
right, as they lend insight into any accretion flow onto a 
binary before merger. A key feature of this flow is the 
presence of spiral density waves in the inner disk cavity. 
Following \3&\ , we highlight the existence of these density 
waves by calculating the surface density fluctuation JS, 
defined by 



(5S = 



(48) 



We also define the torque density, dT/dR, for comparison 
with analytic models and other simulations, 



^ = y" V^T'^.^^.rRdzd^ , 



(49) 



where = (S^)'', which gives = {0,—y,x,0) 
in Cartesian coordinates. Details of the derivation of 
Eq. are given in Appendix [BJ 
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8. Luminosity diagnostics 

In order to study the electromagnetic emission from 
our disk evolutions, we estimate the luminosity due to 
thermal bremmstrahlung and nonthermal synchrotron 
emission using the approximations described in [l6j . For 
synchrotron emission, we assume the presence of a small- 
scale, turbulent B field whose magnitude is approximated 
by setting P = PPm = pB^/{8Tr). We thus assume that 
the magnetic pressure is some fraction 1//3 of its equipar- 
tition value. Simulations of magnetized accretion flows 
have demonstrated that the magnetic fields do not typ- 
ically reach their full equipartition value [Slj. We have 
chosen /3 = 10 to account for this. We also assume that 
the radiation propagates through an optically thin gas 
and we neglect the roles of radiation pressure and radia- 
tive cooling on the hydrodynamic evolution. While an 
accurate estimation of the electromagnetic emission re- 
quires a full solution to the radiative transfer problem, 
this crude method can provide a reasonable estimate of 
the magnitude of the emission under suitable conditions. 

D. Code tests 

Our HRSC general-relativistic hydrodynamic code has 
been thoroughly tested by passing a robust suite of 
tests. These tests include maintaining stable rotating 
stars in stationary equilibrium, reproducing the exact 
Oppenheimer-Snyder solution for collapse to a BH, and 
reproducing analytic solutions for relativistic shocks and 
spherical Bondi accretion onto isolated BHs [tO]. Our 
code has also been used to simulate the collapse of 
very massive, rotating stars to black hole s 18211 : merging 
BHBH binaries [11], BHNS binaries (zl, [sj^ and rela- 
tivistic hydrodynamic matter in the presence of puncture 
black holes [8J| . Recently, our code has been generalized 
to incorporate (optically thick) radiation transport and 
its feedback on fluids in dynamical spacetimes [85|. 

Most of the above tests and simulations were per- 
formed on grids with uniform spacing. In some of the 
simulations, we utilized the multiple-transition fisheye 
transformation (86j so that a uniform computational grid 
spacing corresponds to physical coordinates with spa- 
tially varying resolution. Recently, we have modified our 
code so that we can use the moving-box AMR infrastruc- 
ture provided by Carpet [tH]. To test our new code, we 
have performed shock-tube tests and 3+1 simulations of 
linear gravitational waves, single stationary and boosted 
puncture BHs, puncture BHBH binaries, and rapidly and 
differentially rotating relativistic stars. Our AMR code 
has also been used to perform simulations of BHNS merg- 
ers [i^, binary Bondi and binary Bondi-Hoyle-Lyttleton 
accretion [l^ . 

All of our 3-1-1 AMR code tests were performed as- 
suming equatorial symmetry (i.e., symmetry about the 
z = orbital plane), which we assume in all evolutions 
presented in this paper. We have checked that our AMR 




20 40 60 80 



R / M 



FIG. 1. Surface density profiles at t = (solid line) and 
t ~ tdisk (dotted line), where tdisk is the Keplerian period at 
the radius of maximum pressure. Overlap indicates that disk 
accurately maintains equilibrium configuration over this time 
scale. 

code is able to accurately maintain a stable equilibrium 
disk around a single BH, as demonstrated in Fig. [TJ For 
this test, we use the same disk initial data as run Al in 
Table U except that we set the background metric to be 
that of a single Schwarzchild BH at the origin. As we de- 
scribe in Sec. IIV A[ such a disk is an equilibrium solution 
and is expected to maintain its initial profile. 

We have also checked that the conservations laws in 
Eqs. ((37)) and (|47)) are satisfied in a quasistationary, bi- 
nary spacetime, as described in Sec. IIV CI In Fig. [21 the 
dashed red curve shows the left-hand side of Eq. ([57| . 
with the world tube L chosen to correspond to a sphere 
centered at the orgin with a radius = 25M. For 
comparison, we also plot Mo(t)/Mo,i with the solid red 
curve. The data are from run A2, in which we impose 
a helical Killing vector to solve for the metric as de- 
scribed in Sec. IIII A| while evolving the hydrodynamics 
using Eqs. (HH), (US]), and (ED]). Wc find that Eq. ^ is 
well satisfied, indicating that our code is conserving rest 
mass correctly. Similarly, the dashed black line shows the 
left-hand side of Eq. (H71) . normalized byMg^i, while the 
solid black line shows {E{t) — ilJ{t))/Mo^i for compari- 
son. Again, we see that our code is conserving E — ilJ 
correctly. 



V. RESULTS 

As discussed in Sec. IIIII we separate each of our sim- 
ulations into two phases. We first perform early inspiral 
epoch simulations in which we employ the quasistation- 
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FIG. 2. Plots demonstrating the accurate maintenance of 
conserved quantities. Rest mass Mo{t)/Mo,^ + Jdt {T^p^"^ - 

T^^^)/Mo^i (dashed red line) maintains its initial value accu- 
rately when compared to Mo/Mo,i (solid red). Similarly, the 
conserved quantity {E{t) - nj{t) + J dt (J"^^^ - QTp^) - 
J dt {J-'^'^ — flJ-^^)) /Moj (dashed black line) maintains its 
initial value accurately when compared to {E(t)~nj{t)) /AIo,i 
(solid black line). 



ary CTS metric while keeping the BH separation con- 
stant. This allows the disk to relax to a reliable quasis- 
tationary state. Upon achieving this state, we begin our 
late inspiral and merger epochs simulations in which we 
evolve the metric in full GR, allowing the BHs to inspiral 
and merge. Parameters for each of these disk runs are 
given in Table HI 

Equatorial snapshots from our simulations with F — 
5/3 can be seen in Fig.|3l while meridional snapshots are 
shown in Fig. U) The first two snapshots are from the 
early inspiral epoch calculations, while the second two 
snapshots are from the late inspiral and merger epochs 
calculations. We do not show snapshots for other equa- 
tions of state here, as the accretion flow is qualitatively 
similar. Important results from simulations with other 
equations of state are reported in Table HT] 



A. Early inspiral epoch 

While the formulation outlined in Sec. IIV Al provides 
stable equilibrium disk initial data for a single BH, 
torques from the binary disrupt this equilibrium. Thus, 
it is important to allow the gas to relax to a quasistation- 
ary state before beginning the BH inspiral. We allow this 
relaxation to occur over ~ 5tdisk, where tdisk ~ 1300Af. 
At this time, we find that Mq, L^ym and Lhrem oscillate 



TABLE I. Parameters for BHBH simulations 



Case 


"Epoch 


Orientation 


r 




Al 


early inspiral 


prograde 


5/3 


0.11 


A9 






A 

4/ J 


U.Uo 


A3 






"1.1 


0.03 


A4 




retrograde 


4/3 


0.08 


Bl 


late inspiral 


prograde 


5/3 


0.14 


B2 


and merger 




4/3 


0.11 


B3 






"1.1 


0.06 


B4 




retrograde 


4/3 


0.11 



" Initial binary separation a/M — 10. 

H is the scale height of the disk at ii = Rdiak 
(pressure max), measured at t = for case- A 
runs and a.t t — tmerge for case-B runs. 
" Approximately isothermal. 

around roughly constant values, and there is little evo- 
lution in surface density profiles. Here Lsyn and Lhrem 
refer to synchrotron and bremsstrahlung luminosity, re- 
spectively. 

During the relaxation process, the changes in the mat- 
ter profiles are due to the presence of binary torques. 
These torques cause a disruption of the inner edge of the 
disk, allowing some gas to fall onto the BHs. Infalling gas 
forms spiral waves and shocks which heat the gas near the 
BHs. In the absence of shock heating, the disk would be- 
have adiabatically, and the internal energy density would 
be given by its polytropic value, 

poead = Kpl/{r-l) . (50) 

Thus, shock heating may be measured by computing the 
enhancement in the internal energy density above its adi- 
abatic value and integrating over the disk. We therefore 
compute 

Etnt = / \/^/9oM*e d^x , (51) 

JVa 

and 

Eint,ad = / V-9P0U*<^ad X . (52) 

Here the integral is over Vd, which is the volume between 
r = lOM and the outer boundary of the computational 
domain at I28M. This allows us to compute AEint in 
the bulk of the disk, ignoring the gas near the BHs. In 
Fig. [SI we plot AEint/ Eint,ad vs time during the relax- 
ation of the gas, where AEint = [Ei^t - Eint,ad)- We 
find that AEint /Eint, ad increases monotonically over the 
course of the relaxation process, leveling off to a constant 
value of sa 5 X 10"'' as the disk reaches a quasistation- 
ary state. Because AEint/ Eint, ad <C I, we conclude that 
shock heating does not play a significant role in altering 
the bulk disk profile during this process. 

In order to get a sense of the change in the disk profile 
from the initial data, we plot angle-averaged surface den- 
sity profiles (see Fig. In each case, the torque from 
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FIG. 3. Snapshots of rest-mass density po contours in the orbital plane for cases with F = 5/3. Density contours are plotted at 
po/po,max = ]^Q-2.75+o.5j ^ _ ^ Qy Coutours of highest density are shown with darker shading and are near the BHs. 

Blue arrows denote velocity vectors. The apparent horizon interior is marked by a filled black circle. The top left frame is the 
initial data from the early inspiral epoch, the top right frame is the relaxed, quasistationary disk, which serves as initial data 
for the late inspiral and merger epochs at t « Emerge — 1250Af . The bottom left frame is the data from the late inspiral and 
merger epochs run at i ~ tmerge — 50M, while the bottom left frame is the data from the late inspiral and merger epochs at 

t tmerge ~\~ tdisk' 



the binary has the overall effect of pushing matter out- 
visard. However, the effect of the torques diminishes once 
the disk matter has moved outwrard and a quasistationary 
state is achieved after t ~ Sitdisk- 

In Fig. [21 we plot snapshots of density contours for 
case Al. We see that the disk cavity, which initially ex- 



tends to i? = 15M, becomes partially filled, with a clear 
spiral structure, at small radii R ^ 2a = 20A/ for each 
of the prograde cases. A similar spiral arm structure 
extending inside the disk ca vity has been seen in New- 
tonian simulations (cf. [s^ [sj HI] ) . In the retrograde 
case (A4), this structure is largely absent. This point 
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FIG. 4. Snapshots of rest-mass denstiy po contours in the 
meridional plane for cases with F = 5/3. Density contours 
are plotted at po/po.max = iQ-^ '^s+o-s^ {j = 1,2, ....,6). 
Contours of highest density are near the BHs. The top frame 
is the initial data for the early inspiral epoch; the second 
frame is the relaxed, quasistationary disk that serves as ini- 
tial data for the late inspiral and merger epochs; the third 
frame is the data from the late inspiral and merger epochs at 



t ~ tn 



■ 50M; the bottom frame is the data from the late 



inpiral and merger epochs at f « tmerge + tdisk- 



is further emphasized by comparing the surface density 
fluctuation AS = (S — (S))/ (S) between the prograde 
case A2 and the retrograde case A4 (see Figs. [7] and [5]). 
Following [3^, we compute this quantity in a rotating 
frame in which the binary is stationary, and average over 
several torb- In Fig. [71 we clearly see two strong spiral 
arms emanating from each BH and extending throughout 



t / M 



FIG. 5. Time evolution of lS.Eint / Eint,ad for run A2. 



the cavity region. However, unlike the results of the 2D, 
thin-disk simulations presented in [s^ , we do not see spi- 
ral density waves extending into the bulk of the disk. We 
attribute this to the fact that our 3d disks are thicker, 
which allows waves that are initially propagating in the 
radial direction to be deflected in the vertical direction, 
disrupting the spiral density wave structure. Such effects 
have been demonstrated even for geometrically thin disks 
in which density and temperature are stratified in the z 
direction [Hlgai. 

The structure of the spiral density waves observed near 
Rin is roughly consistent with the theoretical Newtonian 
picture of a wave that is excited by the binary torques 



at the outermost Lindblad resonance at R2 — (3/2)^/'^a 
with orbital resonance VLbin : = 3 : 2 (initially, 
i?2 = 0.87i?i„). We demonstrate this by computing the 
angle-averaged torque density as described in Sec. lIVCTl 
and comparing with the (Newtonian) analytic prediction 
given in Eq. (31) of 3(l|. 



dT 49 o , , . a-* , fR2-B. 
IR - 288-^(^^)"- A^^^ ( ^ 



(53) 



where A2 = 2^-'^{H / R) '^a. Based on numerical data 
from run Al at t > Btdisk, we choose S(i?2) 4 x IO^^'Sq 
and H/R = 0.3, where Eq is the initial maximum surface 
density and H/R is measured at i?2. As we show in 
Fig. [31 we find very good agreement with the analytic 
prediction out to R/a « 2.2, but break down at larger 
radii. This breakdown is not unexpected, as we have 
argued above that spiral density waves do not extend 
into the bulk of the disk as a result of the thickness of 
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FIG. 6. Surface density profiles of tlie F = 5/3 disfc as a 
function of radius. The dashed black curve is the initial 
disk density prole; the solid red one is the surface density 
prole when the disk has reached a quasistationary configura- 
tion after t > 5tdisk, averaged over 2torb- The solid black 
curve is the density profile following the merger, averaged over 
~ 2torb- So is the initial maximum surface density. 

our 3d disks. 

For the retrograde case (A4), we find that the spiral 
density waves are largely absent, as shown in Figs. |S] and 
[TOl This is expected, as the Lindblad resonance does not 
exist when the angular momentum of the disk and the 
binary are antialigned. 

B. Late inspiral and merger epochs 

Having allowed the disk to relax to a quasistationary 
state, we turn to the fully relativistic evolution of met- 
ric and matter fields in order to investigate variations in 
electromagnetic luminosity over the course of the inspiral 
and merger. Our calculations in this epoch apply to the 
decoupling phase of binary inspiral, through merger, but 
before appreciable gas fills the hollow due to viscosity. 
We again consider the prograde cases with F = 4/3 (Bl), 
F = 5/3 (B2), and F = 1.1 (B3) the retrograde case with 
F = 4/3 (B4). In each case, we use the relaxed data from 
the end of the corresponding quasistationary metric run 
as initial data. 

As the binary inspiral proceeds and the separation 
shrinks, the torques due to the binary are diminished. 
As a result, we find that the spiral density waves visible 
when a/M = 10 have largely disappeared by the time of 
merger and remain absent after the merger. This is evi- 
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0/ 27T 

FIG. 7. Time-averaged surface density fluctuation 
(E — (E))/ (E) in the rotating frame in which the binary is at 
rest. The binary point masses are located at R/a = 0.5 and 
</> — (0, tt). Red regions are density maxima and blue regions 
are density minima. Data is from run A2 with F = 4/3. 

dent in the bottom left and right frames of Fig. [3l which 
show snapshots of the density in the equatorial plane 
^ 50M before merger, and ^ Udisk after merger. This 
effect is also illustrated by the evolution of Mq, Li,rem, 
and Lsyn in Fig llll Here, we have computed the lumi- 
nosity assuming a fiducial value of Udisk = 10^^ cm~'^, 
where Udisk is the baryon number density at Rdisk ■ This 
value is consistent with density estimates for a typical 
AGN derived from the Shakura-Sunyaev disk model [oil - 
fosj . albeit in a radiation-dominated, geometrically thin 
regime. However, because there are large variations in 
the gas densities in galactic cores, we provide density scal- 
ings for our results. Because the position of the m = 2 
outermost Lindblad resonance is approximately given by 
i?2 ~ (3/2)^/'^a, we see that as the binary separation is 
reduced, the location of the resonance retreats farther 
inside Rim enabling less matter to be stripped from the 
inner edge of the disk, and reducing Mq. The reduction 
in accretion similarly suppresses the electromagnetic lu- 
minosity generated near the BHs. This effect is exacer- 
bated by the reduction in shock heating due to binary 
torques, as this lowers the temperature of the gas and 
reduces emissivities. Each of these effects is reflected in 
Fig. [TlJ We also show the h+ polarization amplitude 
of the accompanying gravitational wave for comparison. 
Evidently, the decrease in electromagnetic luminosity be- 
ginning at the onset of decoupling is a precursor to the 
late inpiral and merger gravitational radiation. We flnd 
that the choice of EOS can play a significant role in set- 
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FIG. 8. Same as Fig. [7] but for a retrograde disk. 



ting the magnitude of the accretion rate and the lumi- 
nosity. Larger values of F lead to both larger Mq as 
well as higher luminosities. We tabulate the values of lu- 
minosities, accretion rates and characteristic frequencies 
of emission at the onset of the late inspiral and merger 
epochs and just prior to merger in Table HI In addition 
to the increase in the amount of gas near the BHs, larger 
values of F also allow the gas to be shock heated more 
effectively. Because the bremsstrahlung and synchrotron 
emissivities are sensitive to temperature, this also leads 
to an increase in luminosity. Comparing Eq. ()54|) and 
Eq. ([551) below, we see that the temperature dependence 
is much stronger for synchrotron emission. This can ex- 
plain the particularly large differences in synchrotron lu- 
minosity for the different cases reported in Table [III This 
effect also leads the synchrotron luminosity to be domi- 
nated by emission from the heated region near the binary, 
whereas the bremsstrahlung emission is predominantly 
from the bulk of the disk. This dependence accounts 
for the high variability of the synchrotron luminosity in 
comparison to that of the bremsstrahlung emission. 

Because our simulations assume a perfect fluid with no 
dynamical magnetic fields (turbulent fields are assumed 
only to estimate synchrotron emission), there is no viscos- 
ity present to counteract the effect of the binary torques 
in driving matter outward. As a result, we find that 
even after relaxing to a quasistationary disk state in our 
early inspiral epoch calculations, in which the accretion 
rate and luminosity oscillate around fixed values, there 
remains an overall slow outward drift in the bulk of the 
disk. This is evident in Fig. [SI in which the solid red 
curve shows the surface density profile at the beginning 



FIG. 9. Time-averaged torque density dT/dR exerted by 
the binary on the disk after t > 5tdisk- Time averaging was 
carried out over 2torb after the disk has reached a quasis- 
tationary state. Data from run A2 with F = 4/3. The torque 
is plotted in units of lO^'^MaEo. 



of the binary inspiral, while the solid black curve shows 
the surface density profile at t ^ tdisk after the merger. 
We see clear evidence that the bulk of the disk moves out- 
ward, although we suspect that this effect may be altered 
by the inclusion of viscosity. 

In Paper I, we demonstrated that shocks near the BH 
horizons increased in strength throughout the merger as 
the BHs move more supersonically through the surround- 
ing gas. This shock strengthening leads to a tempera- 
ture increase in the inner region, which in turn gives rise 
to an increasing luminosity peaking at the moment of 
merger. Such a temperature increase is largely absent in 
the disklike accretion case treated here, as can be seen 
by comparing temperature contours at the beginning of 
the late inspiral and merger epochs simulation and at 
t ^ tmerge — 50M, as displayed in Fig. [TH As such, we 
do not expect to see significant increase in luminosity 
during the post-decoupling inspiral phase, even at high- 
frequency components of the spectrum. We note that this 
contrasts the conclusions of the Newtonian calculation in 
[l5i] , where a brightening of the precursor light curve be- 
fore merger is found. However, these results do not nec- 
essarily contradict one another, as fll consider geomet- 
rically thin, optically thick disks around non-equal-mass 
BHBH binaries. This is a very different scenario from 
the geometrically thick, optically thin disks surrounding 
equal-mass binaries that we consider in this paper. 

To highlight the role that shock heating plays in our 
simulations, we also plot contours of K/Kq. Here K = 
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FIG. 10. Same as Fig. (9] but for retrograde disk (run A4). 



Pj Pq and Kq is the initial value of K everywhere. The 
quantity K — K(s), where s is the specific gas entropy, 
remains constant in the absence of shocks; shock heating 
yields K/Kq > 1 (see Appendix B of [111). As expected, 
we see that K/Kq increases steeply near the BHs, where 
the binary torque-induced spiral arms are strongest (see 
Fig. I13p . Following the merger, binary torques are no 
longer present and we find that K/Kq is dramatically 
reduced in the cavity region. While we do find that a 
small region near the remnant continues to have K/Kq 
long after the merger (see Fig. [T51) . we note that the gas 
in this region is of very low density and carries a relatively 
insignificant amount of thermal energy. 

We also note that the shocks are confined to the inner 
region and do not propagate into the bulk of the disk. 
While it has been proposed that changes in the potential 
due to mass loss and/or BH kicks following merger can 
give rise to shocks throughout the disk [ll], |93|, we do 
not observe such behavior in our simulations. However, 
this is expected, as it has been noted that a condition for 
shocks to form due to mass loss is that e > H/R, where 
e = {Mi — Mf)/Mi is the fractional mass loss due to grav- 
itational wave emission [ll|. Comparing the fractional 
mass loss for an equal-mass merger, for which e w 0.05, to 
the estimates oi H/R measured at the moment of binary 
merger (see Table U), we find that the above condition is 
never satisfied. The criteria above for shocks to form is 
derived from the condition that the radial velocity must 
exceed the sound speed Cs = (TP/ poh)^^"^ near the inner 
edge of the disk. We have also checked this directly and 
have found that the condition is never met in our simu- 
lations - our disks are too hot, hence too geometrically 



FIG. 11. Time evolution of total BH accretion rate across 
the BH horizons Mo, luminosity L and waveform Dh+ for 
a circumbinary prograde disk with F = 5/3. The initial bi- 
nary separation is a = lOM and the BHs evolve to merger. 
Mq / {pmaxM"^) is the dimensionless accretion rate. Here, 
Pma.M^ = 0.2ni2M|MQyr-i. L = L/[10'"'M|n?2 erg s'^] 
is the total luminosity due to bremsstrahlung (dashed line) 
and synchrotron (solid line) emission. For synchrotron emis- 
sion, we assume /3 = 10. h+ is the -I- polarization of the 
gravitational wave signal as measured by an observer looking 
down the polar axis at a distance D from the binary. BHBH 
merger occurs at t = 0. 



thick, to trigger this effect. 



C. Scaling and detectability 

In quoting values for the accretion rate Mq, we nor- 
malize by the quantity PmaxM'^ — O.2ni2Af|M0yr"^, 
which allows for easy scaling to arbitrary disk density 
and binary mass. It is also possible to derive simple scal- 
ing relations for the luminosities. The dominant region 
of emission differs for bremsstrahlung and synchrotron 
radiation. Because of the stronger dependence of syn- 
chrotron emissivity on temperature (see Appendix B of 
[Tg}). we find that the synchrotron emission originates 
chiefly from the hot gas near the BHs, while the major- 
ity of bremsstrahlung emission originates from the cooler, 
denser gas in the bulk of the disk. In each simulation, 
we find that the temperature is maximum near the BHs, 
typically reaching kTh w lOOMeV at the horizon. In 
the high-temperature limit (kT > irieC^) the synchrotron 



emissivity given in Appendix B of 
perature and density according to 



gsyn oc nlTlP 



11611 scales with tem- 



(54) 
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FIG. 12. Contours of temperature kT (= msP/po) at select 
times for prograde disk with F = 5/3. Contours correspond 
to kT = 940 X W~^+°-^^^ MeV {j = 1,2, ....,6). Frames 
correspond to the beginning of the late inspiral and merger 

epochs phase at t ^ tmerge — 1250M (top), t ~ tmerge — 

50M (middle), and t ~ tmerge + tdisk (bottom). Regions with 
density less than po/po,max < 10~*'^ are left white. Lighter 
shading denotes higher kT. 



FIG. 13. Contours of entropy parameter K (— P/po) at select 
times for prograde disk with F = 5/3. Contours correspond 
to K/Ko = iQO-^+O-^i (j = 1, 2, 6). Frames correspond to 
the beginning of the late inspiral and merger epochs phase at 

t ^ tmerge — 1250M (tOp), t ~ tmerge — 50M (middle), and 

t ~ tmerge + tdiak (bottom). Reglons with dcnslty less than 
po/po,max < 10~*'^ are left white. Darker shading denotes 
higher K. 
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By contrast, the temperature in the bulk of the disk 
where the majority of the bremsstrahlung emission orig- 
inates is nonrelativistic {kTdisk < TigC^) so that the 
bremsstrahlung emissivity scales according to 

Iff « ril,,^T^ll . (55) 

Here and refer to the density and temperature at 
the horizon, and Udisk and Tdisk refer to the density and 
temperature at Rdisk , and (3 = 8nP/ sets the strength 
of the assumed turbulent magnetic field responsible for 
synchrotron emission. As in Paper I, we choose /3 = 10 so 
that the magnetic field is assumed to reach a fraction of 
its equipartition value consistent with MHD simulations 

m- 

Integrating Eqs. ([M]) and ([55|) . we find 

Lsyn ~ J dVqsyn cx nlTlp-^M^ , (56) 

L„./dr,„„„L.r-SflL.W«). (57) 

where the factor of comes from the volume element. 

Because we ignore self-gravity in the disk, our scal- 
ing results apply for arbitrary density. Accordingly, as 
the density at the horizon is varied, its value is sim- 
ply proportional to the maximum density in the disk, 

71^ (X Tldisk ^ Praax- 

Consider the onset of the binary decoupling (late in- 
spiral) phase. The spiral arms through which the gas 
enters the disk cavity are shock heated. The gas in any 
shocked region will be heated to kT ~ niBv'^. Since 
w ^ c near the horizon, shock heating guarantees that 
kTh ^ tubc^ ~ 10^ MeV, independent of the tempera- 
ture in the bulk of the disk (in fact KTh is closer to 100 
MeV) . Scaling for synchrotron luminosity then simplifies 
to 

Lsyn CX ■ (58) 

While the bremsstrahlung luminosity does depend on 
Tdisk, we note that for the fixed values of g, /(_/?.„), and 
Rin/M specified in Eqs. (|M7| . (piS]) . and ((M9)) . the 
enthalpy profile h, and therefore kTdisk — Tnsihdisk ~ 
l)(r — l)/r, is uniquely specified. We therefore regard 
Tdisk as a fixed parameter, so that 

Lff (X pl^axM^ . (59) 

Using these scaling relations, along with the results of 
our simulations, we estimate the average luminosity when 
the binary separation is a = lOM, the adopted onset of 
the late inspiral and merger epochs. Results are given in 
Table HI 

In calculating the luminosity, we have assumed that the 
gas is optically thin. We can verify that this is a good 
assumption by estimating the optical depth. Taking the 
dominant opacity source to be electron scattering, wc find 

Tes ~ uhcttR - 0.2ni2M8 (60) 



where R is the characteristic size of the emission region 
that we have set to Rk, 2M. Thus, we see that our as- 
sumption of an optically thin gas is valid for our canonical 
parameters, although it begins to break down when we 
consider denser disks and/or more massive binaries. 

For bremsstrahlung emission originating at Rdisk, the 
characteristic observed frequency of the emission is given 

by 

hvff kTdisk/ + z) (61) 

for a source at redshift z. We measure the temperature 
at Rdisk for each of our cases, and report the estimated 
characteristic frequencies in Table HT] For canonical pa- 
rameters, bremsstrahlung emission will be predominantly 
in 7 rays. Based on our measured luminosities, we esti- 
mate that the observed flux from this emission will be 
in the range of ~ 10~^^ — 10~^'*n^2-M8^rg cm~^ s~^ for 
a source at 2; — 1. Unfortunately, it is unlikely that 
this emission is strong enough to be detectable. We note 
that the bremsstrahlung emission we measure is actually 
dominated by emission from the bulk of the disk rather 
than the heated gas near the BHs. This makes the emis- 
sion even less likely to be detectable, as it exhibits only 
a small amount of variability. 

In contrast, the synchrotron emission is predominantly 
produced near the BHs. We can estimate the characteris- 
tic frequency of the synchrotron emission by noting that 
Eq. (BIO) of [iBl is maximized when X]\i = 2v/3i^o&^ ~ 
1.09. Here, vq = eB /2TrmeC is the cyclotron frequency 
and 9 = kT/rrieC^. The corresponding observed fre- 
quency is 

1.09 3e/iS { kT V 

hVsyn = — — -; 7 ■ (62) 

1 + z 47rmeC \ meC^ / 

We can use this expression, along with measured values 
of density and temperature in the vicinity of the hori- 
zon, to estimate characteristic values of hvsyn- Values for 
each case at the moment of decoupling and shortly be- 
fore merger are given in Table |TT1 and typically fall in the 
infrared range. We estimate that in each case, the syn- 
chrotron emission should be observable by the proposed 
Wide Field Infrared Survey Telescope (WFIRST) 
and possibly by the Large Synoptic Survey Telescope in- 
strument (LSST) [9^. Our simulations follow the late 
stage of the inspiral in which the binary separation de- 
creases from d = lOM to merger. This corresponds to a 
time scale of 5t ~ lOOMg hrs during which the gradual 
decline in emission should be observed. 



VI. DISCUSSION 

In this paper we have performed a set of fully general- 
relativistic simulations of BHBH binary mergers in a cir- 
cumbinary disk. Our focus has been identifying an ob- 
servable electromagnetic signal that may accompany the 
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TABLE II. Electromagnetic emission at beginning of late inspiral and merger epoch and shortly before merger 



Case 


Time 






^ L f f / 


L syn j L/^Q 


"huff \(1 + zy'^MeV] 


'"fby/L L\_-L 1 12 r^l J 


£31 


i^merge 




n no 
u.uz 


n fit; 




U.o 




B2 






0.003 


0.16 


0.05/3"^ 


0.2 


0.6 


B3 






0.0005 


0.003 


2 X lO'^p-^ 


0.08 


0.01 


B4 






0.01 


0.15 


0.5/3r^ 


0.2 


1.6 


Bl 


tmerge 


- 50M 


0.0003 


0.60 


o.3/?r' 


0.3 


1.7 


B2 






0.0005 


0.14 


o.oosft-^ 


0.2 


1.5 


B3 






0.0001 


0.002 


2.0 X IQ-'^P^^ 


0.08 


1.7 


B4 






0.0005 


0.13 


0.0002/3j:^ 


0.2 


1.4 



"Pma.M^ = O.2ni2M|M0yr-\ ni2 = n/10"cm-^ Mg = M/IO^M©. 
''L46 = 10'"'nf2M|erg s'^ 

■^/^i = /?/10, /? = 8ttP/B\ 



gravitational waves from a black hole merger. Our sim- 
ulations are exploratory only. We have restricted our 
attention to disklike accretion onto equal-mass, nonspin- 
ning BH binaries, although our methods may be extended 
to other binary configurations. We exploit the approx- 
imate helical Killing symmetry to determine the binary 
spacetime for widely separated BHBHs. The disk we 
evolve in this early inspiral spacetime relaxes to near 
quasiequilibrium. Our late inspiral and merger simula- 
tions begin with such a quasistationary disk. We then 
evolve the field as well as the matter. This epoch cor- 
responds to the post disk-binary decoupling phase, ter- 
minating after merger, but before viscosity fills in the 
hollow. 

For each simulation, we have calculated the time- 
varying rest-mass accretion rate, as well as the elec- 
tromagnetic luminosity due to bremsstrahlung and syn- 
chrotron emission. We also derive scaling relations for 
the luminosity, enabling our results to be applied to a 
range of gas parameters and BH masses. 

In each case, we find evidence for a time- varying elec- 
tromagnetic signature accompanying the BHBH binary 
merger. The synchrotron emission is the most easily de- 
tectable component, and we observe a steady decline in 
synchrotron luminosity throughout the post-decoupling 
binary inspiral. This change serves as a characteristic 
precursor of a binary merger, and should be detectable 
by the proposed WFIRST and possibly by the LSST in- 
strument. 

In Paper I, we restricted our attention to Bondi-like 
accetion onto merging binaries. It is instructive to com- 
pare the electromagnetic signatures associated with bi- 
nary Bondi accretion with the signatures from disklike 
accretion discussed in this paper. In the binary Bondi 
case, there is a steady supply of gas accreting onto the 
binary at all stages of the merger. In this case, the evo- 
lution of the luminosity is determined by the strength of 
shock heating near the BHs. As the separation decreases, 
the BHs orbit more rapidly, and the shock-heated tem- 
perature of the gas increases. This increase leads to a 
luminosity that increases throughout the inspiral, then 



drops precipitously following the merger as the shocks 
dissipate. This scenario is quite different from the case 
of disklike accretion, in which binary torques create a hol- 
low region around the binary as well as a small amount 
of matter that leaks into the hollow in the form of spi- 
ral arms. Because the torques decrease throughout the 
inspiral as the BH separation decreases, we find that the 
accretion rate and luminosity decrease steadily over the 
course of our inspiral simulations. 

We suspect, however, that this picture will change with 
the addition of magnetic fields, as the magneto-rotational 
instability will lead to turbulence. We intend to inves- 
tigate this behavior in future calculations, although we 
expect that our results for the late inspiral and merger 
epochs treated here will not be significantly altered. The 
reason is that the time scale for turbulent viscosity to 
fill the hollow with gas for accretion exceeds the inspiral 
time scale following decoupling. 

It is not possible to make a quantitative comparison of 
our results with those of Bode et al [HI due to significant 
differences in our methods. We ernploy disk solutions 
with power-law rotation dependence [T^] and BHBH CTS 
metric data that we allow to relax over a time scale of 
~ 5tdisk ~ 6000Af before we begin our inspiral calcu- 
lations. By contrast. Bode et al employ the constant 
midplane density initial data of [iT] and relax this pro- 
file for a period of ~ 250M. Our calculations also differ 
significantly in that we consider bremmstrahlung radi- 
ation from the entire disk, whereas they consider only 
the cavity region near the BHs. We calculate the syn- 
chrotron emission as well. Nevertheless, we are able to 
see a qualitative agreement in the evolution of the rest- 
mass accretion rate Mq, as a decline in Mq is observed 
throughout the inspiral phase in both calculations. 
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Appendix A: Disk initial data 

Our formulation for an equilibrium stationary disk 
around a single Kerr BH follows closely that of [72, [l^l ■ 

From the conservation of the stress-energy equation, 
we find 



= T^'c-.t 



1 



l3 



"V7 



,/3 



(Al) 



(A2) 



1 



+P,o. + ^(g^^aPohu^u^) . (A3) 

Since we are seeking a solution for a stationary torus in 
Kerr spacetime for a single BH, we can now impose time 
independence, axisymmetry, and no poloidal or radial 
motion; 



dt{...) = d4...) = , 



(A4) 
(A5) 



In Boyer-Lindquist and Kerr-Schild coordinates, these 
constraints imply that Ur = ug = 0. 

We may now simplify Eq. (|A3[) according to 



(A6) 



Here we have also assumed constant entropy, and we have 
introduced the specific angular momentum I = —u^/ut- 
We have also used the fact that u'^Ufj^ = — 1 to show that 

«r' = -(5"-2Z5*^ + ^V*) , (A7) 
and have defined 

n ^ u^/u' = (5*^ - - ig'^) ■ (AS) 

We now follow [t^, [t^ and assume the disk has a power- 
law rotation dependence, whereby takes the form 



17 = ^\-'^ 



where 



<?" - ig' 



(A9) 



(AlO) 



^ 9' 

Combining Eq. (Plj) and Eq. (jAlOl) . we find that 

n = kr , (All) 

where a = q/{q — 2) and k = r/"2/(9-2)_ 



It is now straightforward to show that Eq. (jA6p is sat- 
isfied by 



f{i 



ut{r,0)fil{r,e)) 



(A12) 



where f{l) = |1 - kl''+^\^/^°'+^\ A disk solution is 
uniquely determined for fixed i?i„, Z(i?i„), and q. 

While the solution described above applies to equi- 
librium disks in general Kerr spacetimes, we assume a 
Schwarzschild geometry for the purposes of this study. 

In the Newtonian limit, we know that I = fir^, so 
\^ — l/Q = r^^ whereby 



Thus we can see that asymptotically, 

9 = n= const , 
q = 2 ^ I — const , 
q ~ 1.5 Keplerian 



(A13) 



(A14) 
(A15) 
(A16) 



In this paper, the following parameters are chosen to 
determine the initial disk configuration (set close to Ke- 
plerian) , 



9 = 1.6 , 
1{R„,) = 4.53 , 
= 15M . 



(A17) 
(A18) 
(A19) 



Appendix B: Derivation of the torque density dT/dR 

If we consider a four vector j that is not necessarily 
conserved, then we must generalize Eq. ((27|) to 



dq 



= ~Fh + J"l + 



(Bl) 



If we choose to set H and L to be two concentric cylinders 
of infinite extent, centered around the z axis and of radius 
R and R -t- 5R respectively, we can rewrite Eq. (jBl[) as 

^{q{R + SR) - q{R)) = -T{R) + T{R + SR) 
at 

+ J \/^V^j^i? dR dz d(j) (B2) 
taking the limit 6R — > 0, we therefore find that 

J V^V^^R dz d(j, . 



d I dq 
dR.K'di 



dF 
dR 



We can consider the specific case where j'^ = T^^^cff ^ 
(j)" = {d^Y = (0,-y,a;,0), so that Eq. (|B3| becomes 



(B3) 
and 



d fdJ 
Ir [~dt 



dTt, 



dR 

dT^ 
dR 



J .f^Tt^'V t,(j)vRdzd(j) .(B4) 
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We interpret the first term on the right-hand side of 
Eq. (jB4[) as arising from the net outflow of angular mo- 
mentum carried by matter across the surfaces at R and 
R + SR, while the second term is the torque due to the 
gravitational field. Because we are most interested in 
the torque from the gravitational field of the binary, we 
define 



dT 
dR 



dTt, 



dR 



dF^ 
dR 



^T'"'\/f,(j)^R dz , 



(B5) 



Note that in an axisymmetric spacetime in which (j)^ is a 
Kilhng vector field, T'^'V = T'^''V(^(/)^) = 0, hence 
dT/dR = as expected. 

In order to compute Eq. (IBSP numericaUy, it is conve- 
nient to transform the expression into Cartesian coordi- 
nates. Note that 

= -ry., + T% 

-lyT^''9f.u^. + lxT'^''9^.,y (B6) 



Inserting Eq. (IB6|) into Eq. (IB5|) . we find 
^ = J Rd^ dz^g {-ry, + T%) (B7) 



R 



dzy^ {-yT^^^g^y^x + xT^'^g^^^y) 



Eq. (|B7P is integrated numerically at a number of differ- 
ent radii, so that we can compute profiles of dT/dR. 



Newtonian limit 



We can check that Eq. (jB7p reduces to the correct ex- 
pression in the Newtonian limit. If we let « poi 

\T°i/T^o\ < 1, and \T^^ /T°^\ < 1, we find 



dT _ 1 
dR^ 2 
1 



Rdcj) dzT {~-ygm,x + a;goo,i;) 



- I Rdcl)dz T'^'goo, 



R d(f> dz po$., 



(B8) 



where angled brackets indicate angle averaging, and $ is 
the Newtonian gravitational potential. Eq. (jBSp matches 
the expression given in Eq. (14) of (36j . 
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